home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 22 / Cream of the Crop 22.iso / math / ast53src.zip / PLACALC2.C < prev    next >
C/C++ Source or Header  |  1996-09-29  |  36KB  |  1,002 lines
  1. /*
  2. ** Astrolog (Version 5.30) File: placalc2.c
  3. **
  4. ** IMPORTANT NOTICE: The graphics database and chart display routines
  5. ** used in this program are Copyright (C) 1991-1996 by Walter D. Pullen
  6. ** (Astara@msn.com, http://www.magitech.com/~cruiser1/astrolog.htm).
  7. ** Permission is granted to freely use and distribute these routines
  8. ** provided one doesn't sell, restrict, or profit from them in any way.
  9. ** Modification is allowed provided these notices remain with any
  10. ** altered or edited versions of the program.
  11. **
  12. ** The main planetary calculation routines used in this program have
  13. ** been Copyrighted and the core of this program is basically a
  14. ** conversion to C of the routines created by James Neely as listed in
  15. ** Michael Erlewine's 'Manual of Computer Programming for Astrologers',
  16. ** available from Matrix Software. The copyright gives us permission to
  17. ** use the routines for personal use but not to sell them or profit from
  18. ** them in any way.
  19. **
  20. ** The PostScript code within the core graphics routines are programmed
  21. ** and Copyright (C) 1992-1993 by Brian D. Willoughby
  22. ** (brianw@sounds.wa.com). Conditions are identical to those above.
  23. **
  24. ** The extended accurate ephemeris databases and formulas are from the
  25. ** calculation routines in the program "Placalc" and are programmed and
  26. ** Copyright (C) 1989,1991,1993 by Astrodienst AG and Alois Treindl
  27. ** (alois@azur.ch). The use of that source code is subject to
  28. ** regulations made by Astrodienst Zurich, and the code is not in the
  29. ** public domain. This copyright notice must not be changed or removed
  30. ** by any user of this program.
  31. **
  32. ** Initial programming 8/28,30, 9/10,13,16,20,23, 10/3,6,7, 11/7,10,21/1991.
  33. ** X Window graphics initially programmed 10/23-29/1991.
  34. ** PostScript graphics initially programmed 11/29-30/1992.
  35. ** Last code change made 9/22/1996.
  36. */
  37.  
  38. #include "placalc.h"
  39.  
  40.  
  41. #ifdef PLACALC
  42. /*
  43. ** ---------------------------------------------------------------
  44. ** | Copyright Astrodienst AG and Alois Treindl, 1989,1991,1993  |
  45. ** | The use of this source code is subject to regulations made  |
  46. ** | by Astrodienst Zurich. The code is NOT in the public domain.|
  47. ** |                                                             |
  48. ** | This copyright notice must not be changed or removed        |
  49. ** | by any user of this program.                                |
  50. ** ---------------------------------------------------------------
  51. **
  52. ** Important changes:
  53. ** 11-jun-93 revision 1.12: fixed error which affected Mercury between -2100
  54. ** and -3100 (it jumped wildly).
  55. */
  56.  
  57. #ifdef ASTROLOG
  58. /* Given an object index and a Julian Day time, get its zodiac and        */
  59. /* declination position (planetary longitude and latitude) of the object  */
  60. /* and its velocity and distance from the Earth or Sun. This basically    */
  61. /* just calls the Placalc calculation function to actually do it, but as  */
  62. /* this is the one routine called from Astrolog, this is the one routine  */
  63. /* which has knowledge of and uses both Astrolog and Placalc definitions, */
  64. /* and does things such as translation to Placalc indices and formats.    */
  65.  
  66. bool FPlacalcPlanet(ind, jd, helio, obj, objalt, dir, space)
  67. int ind, helio;
  68. double jd;
  69. real *obj, *objalt, *dir, *space;
  70. {
  71.   int iobj, flag;
  72.   REAL8 jd_ad, rlng, rrad, rlat, rspeed;
  73.  
  74.   if (ind <= oPlu)      /* Convert Astrolog object index to Placalc index. */
  75.     iobj = ind-1;
  76.   else if (ind == oChi)
  77.     iobj = CHIRON;
  78.   else if (FBetween(ind, oCer, oVes))
  79.     iobj = ind - oCer + CERES;
  80.   else if (ind == oNod)
  81.     iobj = us.fTrueNode ? TRUE_NODE : MEAN_NODE;
  82.   else if (ind == oLil)
  83.     iobj = LILITH;
  84.   else
  85.     return fFalse;
  86.  
  87.   jd_ad = jd - JUL_OFFSET;
  88.   flag = helio ? CALC_BIT_SPEED | CALC_BIT_HELIO : CALC_BIT_SPEED;
  89.   jd_ad += deltat(jd_ad);
  90.   if (calc(iobj, jd_ad, flag, &rlng, &rrad, &rlat, &rspeed) == OK) {
  91.     *obj    = rlng;
  92.     *objalt = rlat;
  93.     *dir    = rspeed;
  94.     *space  = rrad;
  95.     return fTrue;
  96.   }
  97.   return fFalse;
  98. }
  99. #endif /* ASTROLOG */
  100.  
  101.  
  102. /***********************************************************
  103. ** $Header$
  104. **
  105. ** definition module for planetary elements
  106. ** and disturbation coefficients
  107. ** version HP-UX C  for new version with stored outer planets
  108. ** 31-jul-88
  109. ** by Alois Treindl
  110. **
  111. ** ---------------------------------------------------------------
  112. ** | Copyright Astrodienst Zurich AG and Alois Treindl, 1989.    |
  113. ** | The use of this source code is subject to regulations made |
  114. ** | by Astrodienst Zurich. The code is NOT in the public domain.|
  115. ** |                                                             |
  116. ** | This copyright notice must not be changed or removed        |
  117. ** | by any user of this program.                                |
  118. ** ---------------------------------------------------------------
  119. **
  120. ***********************************************************/
  121.  
  122. /************************************************************
  123. externally accessible globals, defined as extern in placalc.h
  124. ************************************************************/
  125.  
  126. REAL8 meanekl, ekl, nut;
  127. struct elements el[MARS + 1];
  128.  
  129. /*
  130. ** In the elements degrees were kept as the units for the constants. This
  131. ** requires conversion to radians, when the actual calculations are performed.
  132. ** This approach is not the most efficient, but safer for development.
  133. ** Constant conversion could be done by writing all degree constants with
  134. ** value * DEGTORAD
  135. */
  136.  
  137. #define TIDAL_26  TRUE  /* decide wheter to use new or old lunar tidal
  138. term; a consistent system of delta t must be
  139. used */
  140. #define MOON_TEST_CORR FALSE  /* to include more lunar terms in longitude */
  141.  
  142. REAL8 ekld[4] = {23.452294, -46.845, -.0059, 0.00181};
  143. /* ecliptic with epoch1900, Ekd(0..3) in basic */
  144.  
  145. struct eledata pd[MARS + 1] = {
  146. {/*earth*/  1.00000023, 365.25636042, EPOCH1900,
  147. 99.696678,  .9856473354,  1.089,    0,
  148. 101.220833, 6189.03,  1.63,   0.012,
  149. 0.01675104, -0.00004180,  -0.000000126,
  150. 0,    0,    0,    0,
  151. 0,    0,    0},
  152. /*
  153. ** note 29 June 88 by Alois: G.M.Clemence, Astronomical Journal
  154. ** vol.53,p. 178 (1948) gives a correction to the perihel motion
  155. ** of -4.78" T, giving 6184.25 for the linear Term above. We have
  156. ** not yet applied this correction. It has been used in APAE 22,4
  157. ** on the motion of mars and does make an official impression.
  158. */
  159. {/*moon*/ 0.0025955307, 27.321661,  EPOCH1900,
  160. # if ! TIDAL_26
  161. /*
  162. ** values from Improved Lunar Ephemeris, corresponding to tidal
  163. ** term -22.44"/cy and  consistent with delta t ~ 29.949 T*T
  164. */
  165. 270.4341638,  13.176396526808121, -4.08,  0.0068,
  166. # endif
  167. # if TIDAL_26
  168. /*
  169. ** new values from Morrison 1979, with tidal term -26"/cy as
  170. ** stated in A.E. 1986 onwards, consistent with delta t ~ 44.3 T*T
  171. ** correction: -1.54" + 2.33" T - 1.78" T*T
  172. */
  173. 270.4337361,  13.176396544528099, -5.86,  0.0068,
  174. # endif
  175. 334.329556, 14648522.52,  -37.17,   -0.045,
  176. 0.054900489,  0,    0,
  177. 259.183275, -6962911.23,  7.48 ,    0.008,
  178. 5.145388889,  0,    0},
  179. {/*mercury*/  .3870986, 87.969252,  EPOCH1900,
  180. 178.179078, 4.0923770233, 1.084,    0,
  181. 75.89969722,  5599.76,  1.061,    0,
  182. 0.20561421, .00002046,   -.000000030,
  183. 47.145944444, 4266.75,  .626,   0,
  184. 7.0028805555, 6.699,    -.066},
  185. {/*venus*/  .72333162,  224.700726, EPOCH1900,
  186. 342.767053, 1.6021687039, 1.1148,   0,
  187. 130.16383333, 5068.93,  -3.515,   0,
  188. 0.00682069, -.00004774, .000000091,
  189. 75.7796472223,3239.46,  1.476,    0,
  190. 3.3936305555, 3.621,    .0035},
  191. {/*mars*/ 1.5236914620, 686.9296097,  EPOCH1900,
  192. /* These are the corrected elements by Ross */
  193. 293.74762778, .524071163814,  1.1184,   0,
  194. 334.21820278, 6626.73,  .4675,    -0.0043,
  195. 0.09331290, .000092064, -.000000077,
  196. 48.786441667, 2775.57,  -.005,    -0.0192,
  197. 1.85033333, -2.430,   .0454}
  198. };
  199.  
  200. struct sdat _sd [SDNUM] = {
  201.   114.50,   585.17493,
  202.   109.856,  191.39977,
  203.   148.031,  30.34583,
  204.   284.716,  12.21794,
  205.   114.508,  585.17656,
  206.   -0.56,    359.99213,
  207.   148.03,   30.34743,
  208.   284.72,   12.2196,
  209.   248.07,   1494.726615,
  210.   359.44,   359.993595,
  211.   109.86,   191.402867,
  212.   148.02,   30.348930,
  213.   114.503,  585.173715,
  214.   359.444,  359.989285,
  215.   148.021,  30.344620,
  216.   284.716,  12.21669,
  217.   148.0315, 30.34906264,
  218.   284.7158, 12.22117085,
  219.   220.1695, 4.284931111,
  220.   291.8024, 2.184704167
  221. };
  222.  
  223. REAL8 sa[SDNUM];
  224.  
  225. /*
  226. ** delta long = lampl * COS (lphase - arg) in seconds of arc
  227. ** delta rad  = rampl * COS (rphase - arg) in ninth place of log
  228. ** arg = j * sa (k) + i * ma (this planet)
  229. ** ma = mean anomaly
  230. ** sa = mean anomaly of disturbing planet, where this
  231. ** is taken from the aproximate value in sa[]
  232. ** For the COS (phase - arg) it is good enough to compute
  233. ** with 32 bit reals, because ampl and phase have only
  234. ** four to five significant digits.
  235. ** While saving constant space, it is costing execution time due
  236. ** to float/double conversions.
  237. **
  238. ** In basic, all correction terms for sun, mercury, venus and mars
  239. ** were contained in one array K(0..142,0..6); Nk(N,0) contained
  240. ** the index of the first term of planet N and Nk(N,1) the number
  241. ** of terms for this planet. Here, we use a  0 in the first column
  242. ** kor.j to indicate the end of the table for a planet.
  243. ** K(*) was a basic INTEGER array, therefore the amplitudes and phases
  244. ** had to be expressed as
  245. ** K(i,2) = ampl. of longitude in 0.001 seconds of arc
  246. ** K(i,3) = phase of longitude in 0.01 degrees
  247. ** K(i,4) = ampl. of radius in 9th place of log
  248. ** K(i,5) = phase of radius in 0.01 degrees.
  249. ** Here we have converted the amplitude of long. to seconds of arc
  250. ** and the phases to degrees.
  251. */
  252.  
  253. struct kor ARR earthkor[86+1] = {  /* 11-jul-88 all terms to 0.020" long */
  254.   /* j i  lampl  lphase  rampl  rphase  k */
  255.   -1, 1,  0.013, 243,    28,    335,    8,  /* mercury */
  256.   -1, 3,  0.015, 357,    18,    267,    8,
  257.   -1, 4,  0.023, 326,    5,     239,    8,
  258.   -1, 0,  0.075, 296.6,  94,    205.0,  0,  /* venus */
  259.   -1, 1,  4.838, 299.10, 2359,  209.08, 0,
  260.   -1, 2,  0.074, 207.9,  69,    348.5,  0,
  261.   -1, 3,  0.009, 249,    16,    330,    0,
  262.   -2, 1,  .116,  148.90, 160,   58.40,  0,
  263.   -2, 2,  5.526, 148.31, 6842,  58.32,  0,
  264.   -2, 3,  2.497, 315.94, 869,   226.70, 0,
  265.   -2, 4,  0.044, 311.4,  52,    38.8,   0,
  266.   -3, 2,  0.013, 176,    21,    90,     0,
  267.   -3, 3,  .666,  177.71, 1045,  87.57,  0,
  268.   -3, 4,  1.559, -14.75, 1497,  255.25, 0,
  269.   -3, 5,  1.024, 318.15, 194,   49.50,  0,
  270.   -3, 6,  0.017, 315,    19,    43,     0,
  271.   -4, 4,  .210,  206.20, 376,   116.28, 0,
  272.   -4, 5,  .144,  195.40, 196,   105.20, 0,
  273.   -4, 6,  .152,  -16.20, 94,    254.80, 0,
  274.   -5, 5,  0.084, 235.6,  163,   145.4,  0,
  275.   -5, 6,  0.037, 221.8,  59,    132.2,  0,
  276.   -5, 7,  .123,  195.30, 141,   105.40, 0,
  277.   -5, 8,  .154,  -.40,   26,    270.00, 0,
  278.   -6, 6,  0.038, 264.1,  80,    174.3,  0,
  279.   -6, 7,  0.014, 253,    25,    164,    0,
  280.   -6, 8,  0.01,  230,    14,    135,    0,
  281.   -6, 9,  0.014, 12,     12,    284,    0,
  282.   -7, 7,  0.020, 294,    42,    203.5,  0,
  283.   -7, 8,  0.006, 279,    12,    194,    0,
  284.   -8, 8,  0.011, 322,    24,    234,    0,
  285.   -8, 12, 0.042, 259.2,  44,    169.7,  0,
  286.   -8, 14, 0.032, 48.8,   33,    138.7,  0,
  287.   -9, 9,  0.006, 351,    13,    261,    0,
  288.   1,  -1, .273,  217.70, 150,   127.70, 1,  /* mars */
  289.   1,  0,  0.048, 260.3,  28,    347,    1,
  290.   2,  -3, 0.041, 346,    52,    255.4,  1,
  291.   2,  -2, 2.043, 343.89, 2057,  253.83, 1,
  292.   2,  -1, 1.770, 200.40, 151,   295.00, 1,
  293.   2,  0,  0.028, 148,    31,    234.3,  1,
  294.   3,  -3, .129,  294.20, 168,   203.50, 1,
  295.   3,  -2, .425,  -21.12, 215,   249.00, 1,
  296.   4,  -4, 0.034, 71,     49,    339.7,  1,
  297.   4,  -3, .500,  105.18, 478,   15.17,  1,
  298.   4,  -2, .585,  -25.94, 105,   65.90,  1,
  299.   5,  -4, 0.085, 54.6,   107,   324.6,  1,
  300.   5,  -3, .204,  100.80, 89,    11.00,  1,
  301.   6,  -5, 0.020, 186,    30,    95.7,   1,
  302.   6,  -4, .154,  227.40, 139,   137.30, 1,
  303.   6,  -3, .101,  96.30,  27,    188.00, 1,
  304.   7,  -5, 0.049, 176.5,  60,    86.2,   1,
  305.   7,  -4, .106,  222.70, 38,    132.90, 1,
  306.   8,  -5, 0.052, 348.9,  45,    259.7,  1,
  307.   8,  -4, 0.021, 215.2,  8,     310,    1,
  308.   8,  -6, 0.010, 307,    15,    217,    1,
  309.   9,  -6, 0.028, 298,    34,    208.1,  1,
  310.   9,  -5, 0.062, 346,    17,    257,    1,
  311.   10, -6, 0.019, 111,    15,    23,     1,
  312.   11, -7, 0.017, 59,     20,    330,    1,
  313.   11, -6, 0.044, 105.9,  9,     21,     1,
  314.   13, -8, 0.013, 184,    15,    94,     1,
  315.   13, -7, 0.045, 227.8,  5,     143,    1,
  316.   15, -9, 0.021, 309,    22,    220,    1,
  317.   17, -9, 0.026, 113,    0,     0,      1,
  318.   1,  -2, .163,  198.60, 208,   112.00, 2,  /* jupiter */
  319.   1,  -1, 7.208, 179.53, 7067,  89.55,  2,
  320.   1,  0,  2.600, 263.22, 244,   -21.40, 2,
  321.   1,  1,  0.073, 276.3,  80,    6.5,    2,
  322.   2,  -3, 0.069, 80.8,   103,   350.5,  2,
  323.   2,  -2, 2.731, 87.15,  4026,  -2.89,  2,
  324.   2,  -1, 1.610, 109.49, 1459,  19.47,  2,
  325.   2,  0,  0.073, 252.6,  8,     263,    2,
  326.   3,  -3, .164,  170.50, 281,   81.20,  2,
  327.   3,  -2, .556,  82.65,  803,   -7.44,  2,
  328.   3,  -1, .210,  98.50,  174,   8.60,   2,
  329.   4,  -4, 0.016, 259,    29,    170,    2,
  330.   4,  -3, 0.044, 168.2,  74,    79.9,   2,
  331.   4,  -2, 0.080, 77.7,   113,   347.7,  2,
  332.   4,  -1, 0.023, 93,     17,    3,      2,
  333.   5,  -2, 0.009, 71,     14,    343,    2,
  334.   1,  -2, 0.011, 105,    15,    11,     3,  /* saturn */
  335.   1,  -1, .419,  100.58, 429,   10.60,  3,
  336.   1,  0,  .320,  269.46, 8,     -7.00,  3,
  337.   2,  -2, .108,  290.60, 162,   200.60, 3,
  338.   2,  -1, .112,  293.60, 112,   203.10, 3,
  339.   3,  -2, 0.021, 289,    32,    200.1,  3,
  340.   3,  -1, 0.017, 291,    17,    201,    3,
  341.   ENDMARK
  342. };
  343.  
  344. struct kor ARR mercurykor[24+1] = {
  345.   1,  -1, .711,  35.47,  491,  305.28, 4,
  346.   2,  -3, .552,  161.15, 712,  71.12,  4,
  347.   2,  -2, 2.100, 161.15, 2370, 71.19,  4,
  348.   2,  -1, 3.724, 160.69, 899,  70.49,  4,
  349.   2,  0,  .729,  159.76, 763,  250.00, 4,
  350.   3,  -3, .431,  105.37, 541,  15.53,  4,
  351.   3,  -2, 1.329, 104.78, 1157, 14.84,  4,
  352.   3,  -1, .539,  278.95, 14,   282.00, 4,
  353.   4,  -2, .484,  226.40, 234,  136.02, 4,
  354.   5,  -4, .685,  -10.43, 849,  259.51, 4,
  355.   5,  -3, 2.810, -10.14, 2954, 259.92, 4,
  356.   5,  -2, 7.356, -12.22, 282,  255.43, 4,
  357.   5,  -1, 1.471, -12.30, 1550, 77.75,  4,
  358.   5,  0,  .375,  -12.29, 472,  77.70,  4,
  359.   2,  -1, .443,  218.48, 256,  128.36, 5,
  360.   4,  -2, .374,  151.81, 397,  61.63,  5,
  361.   4,  -1, .808,  145.93, 13,   35.00,  5,
  362.   1,  -1, .697,  181.07, 708,  91.38,  6,
  363.   1,  0,  .574,  236.72, 75,   265.40, 6,
  364.   2,  -2, .938,  36.98,  1185, 306.97, 6,
  365.   2,  -1, 3.275, 37.00,  3268, 306.99, 6,
  366.   2,  0,  .499,  31.91,  371,  126.90, 6,
  367.   3,  -1, .353,  25.84,  347,  295.76, 6,
  368.   2,  -1, .380,  239.87, 0,    0,      7,
  369.   ENDMARK
  370. };
  371.  
  372. struct kor ARR venuskor[22+1] = {
  373.   -1, 2,  .264,   -19.20, 175,  251.10, 8,
  374.   -2, 5,  .361,   167.68, 55,   77.20,  8,
  375.   1,  -1, 4.889,  119.11, 2246, 29.11,  9,
  376.   2,  -2, 11.261, 148.23, 9772, 58.21,  9,
  377.   3,  -3, 7.128,  -2.57,  8271, 267.42, 9,
  378.   3,  -2, 3.446,  135.91, 737,  47.37,  9,
  379.   4,  -4, 1.034,  26.54,  1426, 296.49, 9,
  380.   4,  -3, .677,   165.32, 445,  75.70,  9,
  381.   5,  -5, .330,   56.88,  510,  -33.36, 9,
  382.   5,  -4, 1.575,  193.93, 1572, 104.21, 9,
  383.   5,  -3, 1.439,  138.08, 162,  229.90, 9,
  384.   6,  -6, .143,   84.40,  236,  -5.80,  9,
  385.   6,  -5, .205,   44.20,  256,  314.20, 9,
  386.   6,  -4, .176,   164.30, 70,   75.70,  9,
  387.   8,  -5, .231,   180.00, 25,   75.00,  9,
  388.   3,  -2, .673,   221.62, 717,  131.60, 10,
  389.   3,  -1, 1.208,  237.57, 29,   149.00, 10,
  390.   1,  -1, 2.966,  208.09, 2991, 118.09, 11,
  391.   1,  0,  1.563,  268.31, 91,   -7.60,  11,
  392.   2,  -2, .889,   145.16, 1335, 55.17,  11,
  393.   2,  -1, .480,   171.01, 464,  80.95,  11,
  394.   3,  -2, .169,   144.20, 250,  54.00,  11,
  395.   ENDMARK
  396. };
  397.  
  398. struct kor ARR marskor[62+1] = {
  399.   -1, 1,  .115,   65.84,  684,   156.14, 12,
  400.   -1, 2,  .623,   246.03, 812,   155.77, 12,
  401.   -1, 3,  6.368,  57.60,  556,   -32.06, 12,
  402.   -1, 4,  .588,   57.24,  616,   147.28, 12,
  403.   -2, 5,  .138,   39.18,  157,   309.39, 12,
  404.   -2, 6,  .459,   217.58, 82,    128.10, 12,
  405.   -1, -1, .106,   33.60,  141,   303.45, 13,
  406.   -1, 0,  .873,   34.34,  1112,  304.05, 13,
  407.   -1, 1,  8.559,  35.10,  6947,  304.45, 13,
  408.   -1, 2,  13.966, 20.50,  2875,  113.20, 13,
  409.   -1, 3,  1.487,  22.18,  1619,  112.38, 13,
  410.   -1, 4,  .175,   22.46,  225,   112.15, 13,
  411.   -2, 2,  .150,   18.96,  484,   266.42, 13,
  412.   -2, 3,  7.355,  158.64, 6412,  68.62,  13,
  413.   -2, 4,  4.905,  154.09, 1985,  244.70, 13,
  414.   -2, 5,  .489,   154.39, 543,   244.50, 13,
  415.   -3, 3,  .216,   111.06, 389,   21.10,  13,
  416.   -3, 4,  .355,   110.64, 587,   19.17,  13,
  417.   -3, 5,  2.641,  280.58, 2038,  190.60, 13,
  418.   -3, 6,  .970,   276.06, 587,   6.75,   13,
  419.   -3, 7,  .100,   276.20, 116,   6.40,   13,
  420.   -4, 5,  .152,   232.48, 259,   142.60, 13,
  421.   -4, 6,  .264,   230.47, 387,   139.75, 13,
  422.   -4, 7,  1.156,  41.64,  749,   312.67, 13,
  423.   -4, 8,  .259,   37.92,  205,   128.80, 13,
  424.   -5, 8,  .172,   -8.99,  234,   260.70, 13,
  425.   -5, 9,  .575,   164.48, 308,   74.60,  13,
  426.   -6, 10, .115,   113.70, 145,   23.53,  13,
  427.   -6, 11, .363,   285.69, 144,   196.00, 13,
  428.   -7, 13, .353,   48.83,  85,    319.10, 13,
  429.   -8, 15, 1.553,  170.14, 110,   81.00,  13,
  430.   -8, 16, .148,   170.74, 154,   259.94, 13,
  431.   -9, 17, .193,   293.70, 23,    22.80,  13,
  432.   1,  -3, .382,   46.48,  521,   316.25, 14,
  433.   1,  -2, 3.144,  46.78,  3894,  316.39, 14,
  434.   1,  -1, 25.384, 48.96,  23116, 318.87, 14,
  435.   1,  0,  3.732,  -17.62, 1525,  117.81, 14,
  436.   1,  1,  .474,   -34.60, 531,   59.67,  14,
  437.   2,  -4, .265,   192.88, 396,   103.12, 14,
  438.   2,  -3, 2.108,  192.72, 3042,  102.89, 14,
  439.   2,  -2, 16.035, 191.90, 22144, 101.99, 14,
  440.   2,  -1, 21.869, 188.35, 16624, 98.33,  14,
  441.   2,  0,  1.461,  189.66, 1478,  279.04, 14,
  442.   2,  1,  .167,   191.04, 224,   280.81, 14,
  443.   3,  -4, .206,   167.11, 338,   76.13,  14,
  444.   3,  -3, 1.309,  168.27, 2141,  76.24,  14,
  445.   3,  -2, 2.607,  228.41, 3437,  139.74, 14,
  446.   3,  -1, 3.174,  207.20, 1915,  115.83, 14,
  447.   3,  0,  .232,   207.78, 240,   298.06, 14,
  448.   4,  -4, .178,   127.25, 322,   36.16,  14,
  449.   4,  -3, .241,   200.69, 389,   110.02, 14,
  450.   4,  -2, .330,   267.57, 413,   179.86, 14,
  451.   4,  -1, .416,   221.88, 184,   128.17, 14,
  452.   1,  -2, .155,   -38.20, 191,   231.58, 15,
  453.   1,  -1, 1.351,  -34.10, 1345,  235.85, 15,
  454.   1,  0,  .884,   288.05, 111,   39.90,  15,
  455.   1,  1,  .132,   284.88, 144,   15.67,  15,
  456.   2,  -2, .620,   35.15,  869,   305.30, 15,
  457.   2,  -1, 1.768,  32.50,  1661,  302.51, 15,
  458.   2,  0,  .125,   18.73,  103,   119.90, 15,
  459.   3,  -2, .141,   47.59,  199,   318.06, 15,
  460.   3,  -1, .281,   40.95,  248,   310.75, 15,
  461.   ENDMARK
  462. };
  463.  
  464. struct m45dat m45[NUM_MOON_CORR] = {
  465.   /*    l,     l',  F,  D,      Long,    Lat,         Par),*/
  466.   { 0,   0, 0,   4,    13.902,     14.06,  0.2607},
  467.   { 0,   0, 0,   2,    2369.912,   2373.36,     28.2333},
  468.   { 1,   0, 0,   4,     1.979,      6.98,  0.0433},
  469.   { 1,   0, 0,   2,     191.953,    192.72,  3.0861},
  470.   { 1,   0, 0,   0,   22639.500,  22609.1,     186.5398},
  471.   { 1,   0, 0,  -2,   -4586.465,  -4578.13,     34.3117},
  472.   { 1,   0, 0,  -4, -38.428,    -38.64,  0.6008},
  473.   { 1,   0, 0,  -6,  -0.393,     -1.43,  0.0086},
  474.   { 0,   1, 0,   4,  -0.289,     -1.59, -0.0053},
  475.   { 0,   1, 0,   2, -24.420,    -25.10, -0.3000},
  476.   { 0,   1, 0,   0,    -668.146,   -126.98,     -0.3997},
  477.   { 0,   1, 0,  -2,    -165.145,   -165.06,      1.9178},
  478.   { 0,   1, 0,  -4,  -1.877,     -6.46,  0.0339},
  479.   { 0,   0, 0,   3,   0.403,     -4.01,  0.0023},
  480.   { 0,   0, 0,   1,    -125.154,   -112.79,     -0.9781},
  481.   { 2,   0, 0,   4,   0.213,      1.02,  0.0054},
  482.   { 2,   0, 0,   2,  14.387,     14.78,  0.2833},
  483.   { 2,   0, 0,   0, 769.016,    767.96, 10.1657},
  484.   { 2,   0, 0,  -2,    -211.656,   -152.53,     -0.3039},
  485.   { 2,   0, 0,  -4, -30.773,    -34.07,  0.3722},
  486.   { 2,   0, 0,  -6,  -0.570,     -1.40,  0.0109},
  487.   { 1,   1, 0,   2,  -2.921,    -11.75, -0.0484},
  488.   { 1,   1, 0,   0,    -109.673,   -115.18,     -0.9490},
  489.   { 1,   1, 0,  -2,    -205.962,   -182.36,      1.4437},
  490.   { 1,  1,  0,  -4,  -4.391,     -9.66,  0.0673},
  491.   { 1,  -1, 0,  4,    0.283,      1.53,  0.0060},
  492.   { 1,  -1, 0,  2,   14.577,     31.70,  0.2302},
  493.   { 1,  -1, 0,  0,  147.687,    138.76,  1.1528},
  494.   { 1,  -1, 0,  -2,  28.475,     23.59, -0.2257},
  495.   { 1,  -1, 0,  -4,   0.636,      2.27, -0.0102},
  496.   { 0,  2,  0,  2,   -0.189,     -1.68, -0.0028},
  497.   { 0,  2,  0,  0,   -7.486,     -0.66, -0.0086},
  498.   { 0,  2,  0,  -2,  -8.096,    -16.35,  0.0918},
  499.   { 0,  0,  2,  2,   -5.741,     -0.04, -0.0009},
  500.   { 0,  0,  2,  0,  -411.608,    -0.2,      -0.0124},
  501.   { 0,  0,  2,  -2, -55.173,    -52.14, -0.1052},
  502.   { 0,  0,  2,  -4,   0.025,     -1.67,  0.0031},
  503.   { 1,  0,  0,  1,   -8.466,    -13.51, -0.1093},
  504.   { 1,  0,  0,  -1,  18.609,      3.59,  0.0118},
  505.   { 1,  0,  0,  -3,   3.215,      5.44, -0.0386},
  506.   { 0,  1,  0,  1,   18.023,     17.93,  0.1494},
  507.   { 0,  1,  0,  -1,   0.560,      0.32, -0.0037},
  508.   { 3,  0,  0,  2,    1.060,      2.96,  0.0243},
  509.   { 3,  0,  0,  0,   36.124,     50.64,  0.6215},
  510.   { 3,  0,  0,  -2, -13.193,    -16.40, -0.1187},
  511.   { 3,  0,  0,  -4,  -1.187,     -0.74,  0.0074},
  512.   { 3,  0,  0,  -6,  -0.293,     -0.31,  0.0046},
  513.   { 2,  1,  0,  2,   -0.290,     -1.45, -0.0051},
  514.   { 2,  1,  0,  0,   -7.649,    -10.56, -0.1038},
  515.   { 2,  1,  0,  -2,  -8.627,     -7.59, -0.0192},
  516.   { 2,  1,  0,  -4,  -2.740,     -2.54,  0.0324},
  517.   { 2,  -1, 0,  2,    1.181,      3.32,  0.0213},
  518.   { 2,  -1, 0,  0,    9.703,     11.67,  0.1268},
  519.   { 2,  -1, 0,  -2,  -2.494,     -1.17, -0.0017},
  520.   { 2,  -1, 0,  -4,   0.360,      0.20, -0.0043},
  521.   { 1,  2,  0,  0,   -1.167,     -1.25, -0.0106},
  522.   { 1,  2,  0,  -2,  -7.412,     -6.12,  0.0484},
  523.   { 1,  2,  0,  -4,  -0.311,     -0.65,  0.0044},
  524.   { 1,  -2, 0,  2,    0.757,      1.82,  0.0112},
  525.   { 1,  -2, 0,  0,    2.580,      2.32,  0.0196},
  526.   { 1,  -2, 0,  -2,   2.533,      2.40, -0.0212},
  527.   { 0,  3,  0,  -2,  -0.344,     -0.57,  0.0036},
  528.   { 1,  0,  2,  2,   -0.992,     -0.02,  0},
  529.   { 1,  0,  2,  0,  -45.099,     -0.02, -0.0010},
  530.   { 1,  0,  2,  -2,  -0.179,     -9.52, -0.0833},
  531.   { 1,  0,  -2, 2,   -6.382,     -3.37, -0.0481},
  532.   { 1,  0,  -2, 0,   39.528,     85.13, -0.7136},
  533.   { 1,  0,  -2, -2,   9.366,      0.71,     -0.0112},
  534.   { 0,  1,  2,  0,    0.415,      0.10,    0.0013},
  535.   { 0,  1,  2,  -2,  -2.152,     -2.26, -0.0066},
  536.   { 0,  1,  -2, 2,   -1.440,     -1.30,  0.0014},
  537.   { 0,  1,  -2, -2,   0.384,      0.0,   0.0},
  538.   { 2,  0,  0,  1,   -0.586,     -1.20, -0.0100},
  539.   { 2,  0,  0,  -1,   1.750,      2.01,  0.0155},
  540.   { 2,  0,  0,  -3,   1.225,      0.91,     -0.0088},
  541.   { 1,  1,  0,  1,    1.267,      1.52,  0.0164},
  542.   { 1,  -1, 0,  -1,  -1.089,      0.55,      0},
  543.   { 0,  0,  2,  -1,   0.584,      8.84,  0.0071},
  544.   { 4,  0,  0,  0,    1.938,      3.60,  0.0401},
  545.   { 4,  0,  0,  -2,  -0.952,     -1.58, -0.0130},
  546.   { 3,  1,  0,  0,   -0.551,      0.94, -0.0097},
  547.   { 3,  1,  0,  -2,  -0.482,     -0.57, -0.0045},
  548.   { 3,  -1, 0,  0,    0.681,      0.96,      0.0115},
  549.   { 2,  0,  2,  0,   -3.996,      0,   0.0004},
  550.   { 2,  0,  2,  -2,   0.557,     -0.75,     -0.0090},
  551.   { 2,  0,  -2, 2,   -0.459,     -0.38, -0.0053},
  552.   { 2,  0,  -2, 0,   -1.298,      0.74,      0.0004},
  553.   { 2,  0,  -2, -2,   0.538,      1.14, -0.0141},
  554.   { 1,  1,  -2, -2,   0.426,      0.07,     -0.0006},
  555.   { 1,  -1, 2,  0,   -0.304,      0.03,  0.0003},
  556.   { 1,  -1, -2, 2,   -0.372,     -0.19, -0.0027},
  557.   { 0,  0,  4,  0,    0.418,      0,   0},
  558.   { 2,  -1, 0,  -1,  -0.352,     -0.37,     -0.0028}
  559. };
  560.  
  561. # if MOON_TEST_CORR
  562. /* moon additional correction terms */
  563. struct m5dat {
  564.   REAL8 lng;
  565.   int  i0,i1,i2,i3;
  566. } m5[] = {
  567.   /*    lng,  l,     l',  F,  D,      */
  568.   0.127,  0,  0,  0,  6,
  569.   -0.151, 0,  2,  0,  -4,
  570.   -0.085, 0,  0,  2,  4,
  571.   0.150,  0,  1,  0,  3,
  572.   -0.091, 2,  1,  0,  -6,
  573.   -0.103, 0,  3,  0,  0,
  574.   -0.301, 1,  0,  2,  -4,
  575.   0.202,  1,  0,  -2, -4,
  576.   0.137,  1,  1,  0,  -1,
  577.   0.233,  1,  1,  0,  -3,
  578.   -0.122, 1,  -1, 0,  1,
  579.   -0.276, 1,  -1, 0,  -3,
  580.   0.255,  0,  0,  2,  1,
  581.   0.254,  0,  0,  2,  -3,
  582.   -0.100, 3,  1,  0,  -4,
  583.   -0.183, 3,  -1, 0,  -2,
  584.   -0.297, 2,  2,  0,  -2,
  585.   -0.161, 2,  2,  0,  -4,
  586.   0.197,  2,  -2, 0,  0,
  587.   0.254,  2,  -2, 0,  -2,
  588.   -0.250, 1,  3,  0,  -2,
  589.   -0.123, 2,  0,  2,  2,
  590.   0.173,  2,  0,  -2, -4,
  591.   0.263,  1,  1,  2,  0,
  592.   0.130,  3,  0,  0,  -1,
  593.   0.113,  5,  0,  0,  0,
  594.   0.092,  3,  0,  2,  -2,
  595.   0,  99, 0,  0,  0 /* end mark */
  596. };
  597. # endif /* MOON_TEST_CORR */
  598.  
  599.  
  600. /* solution of the Kepler equation, return rad */
  601. /* t = mean anomaly in degrees                 */
  602. /* ex = excentricity of orbit                  */
  603. /* err = maximum error in degrees              */
  604.  
  605. REAL8 fnu(t, ex, err)
  606. REAL8 t;
  607. REAL8 ex;
  608. REAL8 err;
  609. {
  610.   REAL8 u0, delta;
  611.  
  612.   t *= DEGTORAD;
  613.   u0 = t;
  614.   err *= DEGTORAD;
  615.   delta = 1;
  616.   while (ABS8(delta) >= err) {
  617.     delta = (t + ex * SIN8(u0) - u0) / (1 - ex * COS8(u0));
  618.     u0 += delta;
  619.   }
  620.   return u0;
  621. }
  622.  
  623.  
  624. /* x MOD 360.0, so that x in 0..360 */
  625.  
  626. REAL8 smod8360(x)
  627. REAL8 x;
  628. {
  629.   while (x >= 360.0)
  630.     x -= 360.0;
  631.   while (x < 0.0)
  632.     x += 360.0;
  633.   return x;
  634. }
  635.  
  636.  
  637. /* x MOD 360.0, so that x in 0..360 */
  638.  
  639. REAL8 mod8360(x)
  640. REAL8 x;
  641. {
  642.   if (x >= 0 && x < 360.0)
  643.     return x;
  644.   return x - 360.0 * RFloor(x / 360.0);
  645. }
  646.  
  647.  
  648. /* a - b on a 360 degree circle, result -180..180 */
  649.  
  650. REAL8 diff8360(a, b)
  651. REAL8 a;
  652. REAL8 b;
  653. {
  654.   REAL8 d;
  655.  
  656.   d = a - b;
  657.   if (d >= 180.0)
  658.     return d - 360.0;
  659.   if (d < -180.0)
  660.     return d + 360.0;
  661.   return d;
  662. }
  663.  
  664.  
  665. REAL8 test_near_zero(x)
  666. REAL8 x;
  667. {
  668.   if (ABS8(x) >= NEAR_ZERO)
  669.     return x;
  670.   if (x < 0)
  671.     return -NEAR_ZERO;
  672.   return NEAR_ZERO;
  673. }
  674.  
  675.  
  676. /*
  677. ** p points to memory filled with long values; for
  678. ** each of the values the seqeuence of the four bytes
  679. ** has to be reversed, to translate HP-UX and VAX
  680. ** ordering to MSDOS/Turboc ordering
  681. */
  682.  
  683. void longreorder(p, n)
  684. UCHAR *p;
  685. int n;
  686. {
  687.   int i;
  688.   unsigned char c0, c1, c2, c3;
  689.   static int orderinit = 0;
  690.   unsigned short test;
  691.  
  692.   if (!orderinit) {
  693.     test = 0x0001;
  694.     orderinit = (*(unsigned char *)(&test)) ? 1 : -1;
  695.   }
  696.   if (orderinit < 0)
  697.     return;
  698.   for (i = 0; i < n; i += 4, p += 4) {
  699.     c0 = *p;
  700.     c1 = *(p + 1);
  701.     c2 = *(p + 2);
  702.     c3 = *(p + 3);
  703.     *p = c3;
  704.     *(p + 1) = c2;
  705.     *(p + 2) = c1;
  706.     *(p + 3) = c0;
  707.   }
  708. }
  709.  
  710.  
  711. /*****************************************************
  712. $Header: deltat.c,v 1.10 93/01/27 14:37:06 alois Exp $
  713. deltat.c
  714. deltat(t): returns delta t (in julian days) from universal time t
  715. is included by users
  716. ET = UT +  deltat
  717.  
  718. ---------------------------------------------------------------
  719. | Copyright Astrodienst Zurich AG and Alois Treindl, 1989.    |
  720. | The use of this source code is subject to regulations made  |
  721. | by Astrodienst Zurich. The code is NOT in the public domain.|
  722. |                                                             |
  723. | This copyright notice must not be changed or removed        |
  724. | by any user of this program.                                |
  725. ---------------------------------------------------------------
  726.  
  727. ******************************************************/
  728.  
  729. double deltat(jd_ad)
  730. double jd_ad;
  731. {
  732.   static short int ARR dt[] = { /* in centiseconds */
  733.   /*
  734.   ** dt from 1637 to 2000, as tabulated in A.E.
  735.   ** the values 1620 - 1636 are not taken, as they fit
  736.   ** badly the parabola 25.5 t*t for the next range. The
  737.   ** best crossing point to switch to the parabola is
  738.   ** 1637, where we have fitted the value for continuity
  739.   */
  740.   6780, 6500, 6300,
  741.   6200, 6000, 5800, 5700, 5500,
  742.   5400, 5300, 5100, 5000, 4900,
  743.   4800, 4700, 4600, 4500, 4400,
  744.   4300, 4200, 4100, 4000, 3800, /* 1655 - 59 */
  745.   3700, 3600, 3500, 3400, 3300,
  746.   3200, 3100, 3000, 2800, 2700,
  747.   2600, 2500, 2400, 2300, 2200,
  748.   2100, 2000, 1900, 1800, 1700,
  749.   1600, 1500, 1400, 1400, 1300,
  750.   1200, 1200, 1100, 1100, 1000,
  751.   1000, 1000, 900,  900,  900,
  752.   900,  900,  900,  900,  900,
  753.   900,  900,  900,  900,  900,  /* 1700 - 1704 */
  754.   900,  900,  900,  1000, 1000,
  755.   1000, 1000, 1000, 1000, 1000,
  756.   1000, 1000, 1100, 1100, 1100,
  757.   1100, 1100, 1100, 1100, 1100,
  758.   1100, 1100, 1100, 1100, 1100,
  759.   1100, 1100, 1100, 1100, 1200, /* 1730 - 1734 */
  760.   1200, 1200, 1200, 1200, 1200,
  761.   1200, 1200, 1200, 1200, 1300,
  762.   1300, 1300, 1300, 1300, 1300,
  763.   1300, 1400, 1400, 1400, 1400,
  764.   1400, 1400, 1400, 1500, 1500,
  765.   1500, 1500, 1500, 1500, 1500, /* 1760 - 1764 */
  766.   1600, 1600, 1600, 1600, 1600,
  767.   1600, 1600, 1600, 1600, 1600,
  768.   1700, 1700, 1700, 1700, 1700,
  769.   1700, 1700, 1700, 1700, 1700,
  770.   1700, 1700, 1700, 1700, 1700,
  771.   1700, 1700, 1600, 1600, 1600, /* 1790 - 1794 */
  772.   1600, 1500, 1500, 1400, 1400,
  773.   1370, 1340, 1310, 1290, 1270, /* 1800 - 1804 */
  774.   1260, 1250, 1250, 1250, 1250,
  775.   1250, 1250, 1250, 1250, 1250,
  776.   1250, 1250, 1240, 1230, 1220,
  777.   1200, 1170, 1140, 1110, 1060,
  778.   1020, 960,  910,  860,  800,
  779.   750,  700,  660,  630,  600,  /* 1830 - 1834 */
  780.   580,  570,  560,  560,  560,
  781.   570,  580,  590,  610,  620,
  782.   630,  650,  660,  680,  690,
  783.   710,  720,  730,  740,  750,
  784.   760,  770,  770,  780,  780,
  785.   788,  782,  754,  697,  640,  /* 1860 - 1864 */
  786.   602,  541,  410,  292,  182,
  787.   161,  10, -102, -128, -269,
  788.   -324, -364, -454, -471, -511,
  789.   -540, -542, -520, -546, -546,
  790.   -579, -563, -564, -580, -566,
  791.   -587, -601, -619, -664, -644, /* 1890 - 1894 */
  792.   -647, -609, -576, -466, -374,
  793.   -272, -154, -2, 124,  264,
  794.   386,  537,  614,  775,  913,
  795.   1046, 1153, 1336, 1465, 1601,
  796.   1720, 1824, 1906, 2025, 2095,
  797.   2116, 2225, 2241, 2303, 2349, /* 1920 - 1924 */
  798.   2362, 2386, 2449, 2434, 2408,
  799.   2402, 2400, 2387, 2395, 2386,
  800.   2393, 2373, 2392, 2396, 2402,
  801.   2433, 2483, 2530, 2570, 2624,
  802.   2677, 2728, 2778, 2825, 2871,
  803.   2915, 2957, 2997, 3036, 3072, /* 1950 - 1954 */
  804.   3107, 3135, 3168, 3218, 3268,
  805.   3315, 3359, 3400, 3447, 3503,
  806.   3573, 3654, 3743, 3829, 3920,
  807.   4018, 4117, 4223, 4337, 4449,
  808.   4548, 4646, 4752, 4853, 4959,
  809.   5054, 5138, 5217, 5296, 5379, /* 1980 - 1984 */
  810.   5434, 5487, 5532, 5582, 5630, /* 1985 - 89 from AE 1991 */
  811.   5686, 5757, 5900, 5900, 6000, /* AE 1993 and extrapol */
  812.   6050, 6100, 6150, 6200, 6250, /* 1995 - 1999 */
  813.   6300};          /* 2000 */
  814.   double yr, cy, delta;
  815.   long iyr, i;
  816.   yr = (jd_ad + 18262) / 365.25 + 100.0;    /*  year  relative 1800 */
  817.   cy = yr / 100;
  818.   iyr =  (long) (RFloor(yr) + 1800);   /* truncated to integer, rel 0 */
  819. #if TIDAL_26    /* Stephenson formula only when 26" tidal
  820.   term in lunar motion */
  821.   if (iyr >= 1637  && iyr < 2000) {
  822.   i = iyr - 1637;
  823.   delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - RFloor(yr)) * 0.01;
  824.   } else if (iyr >= 2000) { /* parabola, fitted at value[2000] */
  825.   delta = 25.5 * cy * cy  - 25.5 * 4 + 63.00;
  826.   } else if (iyr >= 948) {  /* from 948 - 1637 use parabola */
  827.   delta = 25.5 * cy * cy;
  828.   } else {  /* before 984 use other parabola */
  829.   delta = 1361.7  + 320 * cy + 44.3 * cy * cy;  /* fits at 948 */
  830.   }
  831. #else    /* use Clemence value + 5 sec before 1690, new dt afterwards */
  832.   cy -= 1;  /* epoch 1900 */
  833.   if (iyr >= 1690  && iyr < 2000) {
  834.   i = iyr - 1637;
  835.   delta = dt[i] * 0.01 + (dt[i+1] - dt[i]) * (yr - RFloor(yr)) * 0.01;
  836.   } else if (iyr >= 2000) { /* parabola, fitted at value[2000] */
  837.   delta = 29.949 * cy * cy  - 29.949 * 4 + 63.0;
  838.   } else {
  839.   delta = 5 + 24.349 + 72.3165 * cy + 29.949 * cy * cy; /* fits at 1690 */
  840.   }
  841. #endif
  842.   return delta / 86400.0;
  843. }
  844.  
  845.  
  846. /*******************************************
  847. $Header: d2l.c,v 1.9 91/11/16 16:24:20 alois Exp $
  848.  
  849. double to long with rounding, no overflow check
  850. *************************************/
  851. long d2l(x)
  852. double x;
  853. {
  854.   if (x >=0)
  855.     return ((long) (x + 0.5));
  856.   else
  857.     return (-(long) (0.5 - x));
  858. }
  859.  
  860.  
  861. /*
  862. * $Header$
  863. *
  864. * A collection of useful functions for centisec calculations.
  865.  
  866. ---------------------------------------------------------------
  867. | Copyright Astrodienst Zurich AG and Alois Treindl, 1991.    |
  868. | The use of this source code is subject to regulations made  |
  869. | by Astrodienst Zurich. The code is NOT in the public domain.|
  870. |                                                             |
  871. | This copyright notice must not be changed or removed        |
  872. | by any user of this program.                                |
  873. ---------------------------------------------------------------
  874. *******************************************************/
  875.  
  876. double degnorm(p)
  877. double p;
  878. {
  879.   if (p < 0)
  880.     do {
  881.       p += 360.0;
  882.     } while (p < 0);
  883.   else if (p >= 360.0)
  884.     do {
  885.       p -= 360.0;
  886.     } while (p >= 360.0);
  887.   return (p);
  888. }
  889.  
  890.  
  891. /*********************************************************
  892. $Header: julday.c,v 1.9 91/11/16 16:25:06 alois Exp $
  893. *********************************************************/
  894.  
  895. /*
  896. ** This function returns the absolute Julian day number (JD)
  897. ** for a given calendar date.
  898. ** The aruguments are a calendar date: day, month, year as integers,
  899. ** hour as double with decimal fraction.
  900. ** If gregflag = 1, Gregorian calendar is assumed, gregflag = 0
  901. ** Julian calendar is assumed.
  902. **
  903. ** The Julian day number is system of numbering all days continously
  904. ** within the time range of known human history. It should be familiar
  905. ** for every astrological or astronomical programmer. The time variable
  906. ** in astronomical theories is usually expressed in Julian days or
  907. ** Julian centuries (36525 days per century) relative to some start day;
  908. ** the start day is called 'the epoch'.
  909. ** The Julian day number is a double representing the number of
  910. ** days since JD = 0.0 on 1 Jan -4712, 12:00 noon.
  911. ** Midnight has always a JD with fraction .5, because traditionally
  912. ** the astronomical day started at noon.
  913. **
  914. ** NOTE: The Julian day number is named after the monk Julianus. It must
  915. ** not be confused with the Julian calendar system, which is named after
  916. ** Julius Cesar, the Roman politician who introduced this calendar.
  917. **
  918. ** Original author: Marc Pottenger, Los Angeles.
  919. ** with bug fix for year < -4711   15-aug-88 by Alois Treindl
  920. **
  921. ** References: Oliver Montenbruck, Grundlagen der Ephemeridenrechnung,
  922. **             Verlag Sterne und Weltraum (1987), p.49 ff
  923. **
  924. ** related functions: revjul() reverse Julian day number: compute the
  925. **             calendar date from a given JD
  926. */
  927.  
  928. double julday(month, day, year, hour, gregflag)
  929. int month;
  930. int day;
  931. int year;
  932. double hour;
  933. int gregflag;
  934. {
  935.   double jd, u, u0, u1, u2;
  936.  
  937.   u = year;
  938.   if (month < 3)
  939.     u -=1;
  940.   u0 = u + 4712.0;
  941.   u1 = month + 1.0;
  942.   if (u1 < 4)
  943.     u1 += 12.0;
  944.   jd = RFloor(u0*365.25)
  945.     + RFloor(30.6*u1+0.000001)
  946.     + day + hour/24.0 - 63.5;
  947.   if (gregflag) {
  948.     u2 = RFloor(ABS8(u) / 100) - RFloor(ABS8(u) / 400);
  949.     if (u < 0.0)
  950.       u2 = -u2;
  951.     jd = jd - u2 + 2;
  952.     if ((u < 0.0) && (u/100 == RFloor(u/100)) && (u/400 != RFloor(u/400)))
  953.       jd -= 1;
  954.   }
  955.   return jd;
  956. }
  957.  
  958.  
  959. /*********************************************************
  960. $Header: revjul.c,v 1.9 91/11/16 16:25:37 alois Exp $
  961. *********************************************************/
  962.  
  963. /*
  964. ** revjul() is the inverse function to julday(), see the description there.
  965. ** Arguments are julian day number, calendar flag (0=julian, 1=gregorian)
  966. ** return values are the calendar day, month, year and the hour of
  967. ** the day with decimal fraction (0 .. 23.999999).
  968. **
  969. ** Original author Mark Pottenger, Los Angeles.
  970. ** with bug fix for year < -4711 16-aug-88 Alois Treindl
  971. */
  972.  
  973. void revjul(jd, gregflag, jmon, jday, jyear, jut)
  974. double jd;
  975. int gregflag;
  976. int *jmon;
  977. int *jday;
  978. int *jyear;
  979. double *jut;
  980. {
  981.   double u0, u1, u2, u3, u4;
  982.  
  983.   u0 = jd + 32082.5;
  984.   if (gregflag) {
  985.     u1 = u0 + RFloor(u0/36525.0) - RFloor(u0/146100.0) - 38.0;
  986.     if (jd >= 1830691.5) u1 +=1;
  987.       u0 = u0 + RFloor(u1/36525.0) - RFloor(u1/146100.0) - 38.0;
  988.   }
  989.   u2 = RFloor(u0 + 123.0);
  990.   u3 = RFloor((u2 - 122.2) / 365.25);
  991.   u4 = RFloor((u2 - RFloor(365.25 * u3)) / 30.6001);
  992.   *jmon = (int)(u4-1.0);
  993.   if (*jmon > 12)
  994.     *jmon -= 12;
  995.   *jday = (int)(u2 - RFloor(365.25 * u3) - RFloor(30.6001 * u4));
  996.   *jyear = (int)(u3 + RFloor((u4 - 1.9999) / 12.0) - 4800.0);
  997.   *jut = (jd - RFloor(jd + 0.5) + 0.5) * 24.0;
  998. }
  999. #endif /* PLACALC */
  1000.  
  1001. /* placalc2.c */
  1002.